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Methods and Devices for CT Reconstruction Using a Grangeat Approach 

CROSS REFERENCE TO RELATED APPLICATION 

This application claims priority to U.S. Provisional Application No. 60/447,580 
filed February 14, 2003, and hereby incorporated by reference. 

STATEMENT REGARDING FEDERALLY 
SPONSORED RESEARCH OR DEVELOPMENT 

This invention was made with government support under NIH Grants ROl 

DC03590 and EB001685. The government may have certain rights in the invention. 

BACKGROUND 

This invention relates to methods and devices for volumetric computed 
tomograph (CT) reconstruction and more specifically to methods and devices for 
Grangeat-type half-scan cone beam volumetric CT reconstruction for short and long 
objects in the helical and circular scanning cases. 

In at least one known CT imaging system configuration, an x-ray source 
projects a fan-shaped beam which is collimated to lie within an X-Y plane of a 
Cartesian coordinate system, generally referred to as the "imaging plane". The x-ray 
beam passes through the object being imaged, such as a patient. The beam, after being 
attenuated by the object, impinges upon an array of radiation detectors. The intensity 
of the attenuated beam radiation received at the detector array is dependent upon the 
attenuation of the x-ray beam by the object. Each detector element of the array 
produces a separate electrical signal that is a measurement of the beam attenuation at 
the detector location. The attenuation measurements from all the detectors are acquired 
separately to produce a transmission profile. 

In other known CT systems, the x-ray source and the detector array are rotated 
with a gantry within the imaging plane and around the object to be imaged so that the 
angle at which the x-ray beam intersects the object constantly changes. A group of 
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X ray attenuation measurements, i.e., projection data, from the detector array at one 
gantry angle is referred to as a "view". A "scan" of the object comprises a set of views 
made at different gantry angles, or view angles, during one revolution of the x-ray 
source and detector. In an axial scan, the projection data is processed to construct an 
5 image that corresponds to a two dimensional slice taken through the object. 

One method for reconstructing an image from a set of projection data is referred 
to in the art as the filtered back-projection technique. This process converts the 
attenuation measurements from a scan into integers called "CT numbers" or 
"Hounsfield units", which can be used to control the brightness of a corresponding 

10 pixel on a cathode ray tube display. 

The two-dimensional methods discussed above can reconstruct a slice of the 
measured object. If a volume segment needs to be reconstructed, the complete 
procedure can be performed slice-by-slice with a small movement of the object or of 
the source-detector system between each slice. 

15 A more efficient acquisition setup for volumetric CT uses a two-dimensional 

detector. The rays then form a cone with its base on the detector and its apex on the 
source. An x-ray source naturally produces a cone of rays, so cone-beam acquisition 
not only increases the scanning speed, but also makes better use of the emitted rays 
otherwise wasted by collimation. 

20 Modem CT scanners are rapidly moving from fan-beam towards cone-beam 

geometry. Current micro-CT scanners are already in cone-beam geometry. Half-scan 
CT algorithms are advantageous in terms of temporal resolution and are widely used in 
fan-beam and cone-beam geometry. 

A helical source trajectory is natural for volume scanning of long objects. A 

25 continuously translated object and a rotating source-detector system yield a helical 
source trajectory around the object. Helical scanning has been used for many years 
with one-dimensional detectors and has now been extended for use with multi-row 
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detectors with potential applications for two-dimensional detectors in the medical 
imaging field. 

Sixteen-slice helical computerized tomography scanners are conmaercially 
available. An efficient way to acquire volumetric patient data is to use a helical source 
trajectory together with a multi-row detector. 

With the increasing number of detector rows, the cone-beam CT systems are 
expected to be available in the near future. New applications made possible by these 
new fast volumetric imaging technologies include, but are not limited to, cardiac and 
lung examinations, CT angiography, and interventional procedures. In those 
applications, high temporal resolution is one of the most important requirements. 

Half-scan techniques have been developed to improve the temporal resolution 
for axial and spiral CT images. Various half-scan cone-beam reconstruction methods 
are known in the art. 

The first practical algorithm for three-dimensional reconstruction from cone 
beam projections acquired from a circular source trajectory was the Feldkamp method. 
Various Feldkamp-type algorithms were developed for half-scan cone-beam 
reconstructions from half-scan data collected from circular and helical loci. This 
method has certain limitations, most notably off-mid-plane artifacts that occur because 
the approximation error becomes larger as it goes away from the mid-plane. 

It is therefore desirable to provide a method and apparatus for a half-scan cone- 
beam, three dimensional reconstruction of computed tomographic images that performs 
appropriate data filling using the Grangeat approach and suppresses the off-mid-plane 
artifacts associated with the Feldkamp-type algorithms. 

SUMMARY 

The present invention is directed to methods and devices for volumetric CT 
reconstruction. According to an exemplary embodiment, the device may include a 
system processor that supports the desired functionality as described in detail below 
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and a system data store (SDS) that stores data associated with this functionality, such as 
Grangeat-type half-scan cone beam volumetric CT reconstruction data. 

The system processor is in communication with the SDS. The SDS may include 
multiple physical and/or logical data stores for storing the various types of information 
5 used. Data storage and retrieval functionality can be provided by either the system 
processor or one or more data storage processors associated with the SDS. The system 
processor may communicate with the SDS via any suitable communication channel(s). 
The system processor may include one or more processing elements that are adapted or 
programmed to support the desired image reconstruction and/or other functionality. 

10 Currently, cone-beam CT and Micro-CT scanners are under rapid development 

for major biomedical applications. Half-scan cone-beam image reconstruction 
algorithms assume only part of a scanning turn, and are advantageous in terms of 
temporal resolution and image artifacts. While half-scan cone-beam algorithms known 
in the art are in the Feldkamp framework, according to exemplary embodiments a half- 

15 scan algorithm in the Grangeat framework is used for circular and helical trajectories. 
According to exemplary embodiments, a method of image reconstruction 
includes a variety of steps that may, in certain embodiments, be executed by the 
environment summarized above and more fully described below or be stored as 
computer executable instructions in and/or on any suitable combination of computer- 

20 readable media. 

According to one embodiment, a half-scan algorithm is used for image 
reconstruction for a circular scanning case. The half-scan spans 180 degrees plus two 
cone angles that provide sufficient data for reconstruction of the mid-plane defined by 
the source trajectory. A generalized half-scan, such as a so-called extended half scan, 

25 is also possible. Smooth half-scan weighting functions are used for suppression of data 
inconsistencies in some embodiments. Numerical simulation results are reported for 
verification of formulas and programs of the present invention. This Grangeat-type 
half-scan algorithm produces excellent image quahty, without off-mid-plane artifacts 
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associated with Feldkamp-type half-scan algorithms. The Grangeat-type half-scan 
algorithm can be used for quantitative and dynamic biomedical CT and micro-CT 
applications. 

In another embodiment, image reconstruction may be implemented in a helical 
case without data truncation. According to this embodiment, Grangeat's formula is 
modified for utilization and estimation of Radon data. Specifically, each characteristic 
point in the Radon space can be categorized into singly, doubly, triply sampled, and 
shadow regions. A smooth weighting strategy can be used to compensate for data 
redundancy and inconsistency in some embodiments. 

In the helical half-scan case, the projected trajectories and transition points on 
meridian planes are introduced to guide the design of weighting functions. Then, the 
shadow region is recovered via linear interpolation after smooth weighting. The 
Shepp-Logan phantom and other phantoms can be used to verify the correctness of the 
formulation. The Grangeat-type helical half-scan algorithm of the present invention is 
not only valuable for quantitative and/or dynamic biomedical applications of CT and 
micro-CT, but it can also serve as a basis for solving the long object problem in the 
Grangeat framework. 

Additional advantages of the invention will be set forth in part in the description 
which follows, and in part will be obvious from the description, or may be learned by 
practice of the invention. It is to be understood that both the foregoing general 
description and the following detailed description and the description in the 
attachments hereto are exemplary and explanatory only and are not restrictive of the 
invention. 

BRIEF DESCRIPTION OF THE DRAWINGS 

The accompanying drawings, which are incorporated in and constitute a part of 
this specification, illustrate embodiments of the invention and together with the 
description, serve to explain the principles of the invention. 
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FIG. 1 illustrates an exemplary device and system for image reconstruction of 
an object; 

FIG. 2 illustrates exemplary cone beam geometry and plane integration which is 
equivalent to a 3 Radon transform according to exemplary embodiments; 
5 FIGS. 3(a)-(c) illustrate exemplary meridian and detector planes for a circular 

scan. FIG. 3(a) shows a relationship between meridian and detector planes, FIG. 3(b) 
illustrates a meridian plane, and FIG. 3(c) illustrates a detector plane; 

FIG. 4 illustrates exemplary half-scan geometry as it relates to a first 
embodiment; 

10 FIGS. 5(a)-(d) illustrate a number of exemplary measured detector planes and 

critical angles in the four angular intervals II, 12, 13, and 14, respectively; 
FIG. 6 illustrates exemplary plots of the critical angle functions; 
FIGS. 7(a) and 7(b) illustrate exemplary shadow zone geometry for a mid-plane 
view and a meridian plane view, respectively; 
15 FIGS. 8(a)-8(d) illustrate exemplary maps of singly and doubly sampled zones, 

as well as shadow zones on different meridian planes. The angles represented in 
FIGS. 8(a)-8(d) are 7, 40, 103, and 170 degrees, respectively; 

FIGS. 9(a) and 9(b) illustrate exemplary smooth weighting functions for the 
four maps shown in FIGS. 8(a)-8(d); 
20 FIGS. 10(a) and 10(b) illustrate an exemplary first derivative Radon data of the 

3D Shepp-Logan phantom after data filling with zero padding and linear interpolation, 
respectively; 

FIGS. 1 1(a) and 11(b) illustrate exemplary Grangeat-type half -scan 
reconstruction of the sphere phantom in the plane y=0 with zero-padding and linear 
25 interpolation, respectively; 

FIG. 1 1(c) illustrates exemplary comparative profiles from the identical 
positions marked by the dashed line in 1 1(a); 
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FIGS. 12(a) and 12(b) illustrate exemplary Grangeat-type half-scan 
reconstruction of the 3D Shepp-Logan phantom with zero-padding and linear 
interpolation, respectively; 

FIGS. 13(a) and 13(b) illustrate exemplary Grangeat-type reconstruction of the 
5 3D Shepp-Logan phantom with full-scan reconstruction and half-scan reconstruction, 
respectively; 

FIG. 13(c) illustrates exemplary comparative profiles from the identical 
positions marked by the dashed line in FIG. 13(a); 

FIGS. 14(a)-(c) illustrate exemplary Feldkamp-type half-scan and Grangeat- 
10 type half-scan reconstructions with the original phantom, Feldkamp-type half-scan 
reconstruction, and Grangeat-type half-scan reconstruction, respectively; 

FIG. 14(d) illustrates exemplary comparative profiles from the identical 
positions marked by the dashed line in FIG. 14(a); 

FIG. 15 depicts another illustration of exemplary Grangeat cone beam geometry 
15 according to exemplary embodiments; 

FIGS. 16(a)-(c) depict another illustration of exemplary meridian and detector 
planes with a helical scanning locus explicitly drawn. FIG. 16(a) shows a relationship 
between meridian and detector planes, FIG. 16(b) illustrates a meridian plane, and 
FIG. 16(c) illustrates a detector plane, all of which are in the context of a helical 
20 scanning locus; 

FIG. 17 illustrates exemplary helical half-scan geometry as it relates to the 
second embodiment; 

FIGS. 18(a)-(d) illustrate exemplary classification of the Radon space in the 
shadow region, singly sampled region, doubly sampled region, and triply sampled 
25 region, respectively; 

FIGS. 19(a)-19(i) illustrate exemplary projected trajectories on meridian planes 
of ^= 90^ 110°, 130°, 150°, 170°, 190°, 210°, 230°, and 250° respectively; 
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FIGS. 20(a)-(d) illustrate exemplary transition points for design of the 
weighting functions on meridian planes are shown, and in FIGS. 20(b) - (d), case i, 
case ii, and case viii are shown, respectively; 

FIGS. 21(a)-(d), FIGS. 21(eHh), and FIGS. 21(i)-(l) illustrate exemplary region 
map and weighting functions for meridian planes at ^ = 90% ^ = 170" , and (p - 250** , 
respectively. In FIGS. 21(a), (e), (i), the region map is the brightest is for triply 
sampled zones and the darkest shadow zone. In FIGS. 21(b), (f), 0)> the region map is 
the brightest for the first weighting distribution. In FIGS. 21(c), (g), (k), the region 
map is the brightest for the second weighting distribution, and in FIGS. 21(d), (h), (1), 
the region map is the brightest for the third weighting distribution; 

FIGS. 22(a) and (b) illustrate an exemplary first derivative radon data of the 3D 
Shepp-Logan phantom after data filling with zero padding and linear interpolation, 
respectively; 

FIGS. 23(a)-(g) illustrate an exemplary first derivative Radon data of the 3D 
Shepp-Logan phantom at <p^9^ for (a) group 1, (b) group 2, and (c) group3, (d)(e)(f) 
weighting functions for each group and (g) combined data; 

FIGS. 24(a) and (b) illustrate exemplary reconstructed images of the 3D Shepp- 
Logan phantom for helical half-scan Feldkamp and Grangeat-type helical half-scan 
reconstruction with zero-padding and with linear interpolation, respectively; and 

FIGS. 25(a) and (b) illustrate exemplary long object image reconstruction. 

DETAILED DESCRIPTION 

One or more exemplary embodiments are now described in detail hereinbelow 
and illustrated in the attached drawings. Referring to the drawings, like numbers 
indicate like parts throughout the views. As used in the description herein, the meaning 
of "a," "an," and "the** includes plural reference unless the context clearly dictates 
otherwise. Also, as used in the description herein, the meaning of "in** includes "in" 
and "on" unless the context clearly dictates otherwise. Finally, as used in the 
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description herein, the meanings of "and" and "or" include both the conjunctive and 
disjunctive and may be used interchangeably unless the context clearly dictates 
otherwise. 

Ranges may be expressed herein as from "about" one particular value, and/or to 
5 "about" another particular value. When such a range is expressed, another embodiment 
includes from the one particular value and/or to the other particular value. Similarly, 
when values are expressed as approximations, by use of the antecedent "about," it will 
be understood that the particular value forms another embodiment. It will be further 
understood that the endpoints of each of the ranges are significant both in relation to the 

10 other endpoint, and independentiy of the other endpoint. 

According to an exemplary embodiment, the invention can be implemented in 
conjunction with an x-ray tomographic system, such as that depicted schematically in 
FIG. 1. While an exemplary system using x-rays is described below, one skilled in the 
art will appreciate that the cone-beam algorithms could be used, with or without 

15 modification as necessary, in conjunction with other imaging modalities that involve 
cone-beam geometry, including but not limited to PET and SPECT scanners. 

As shown in FIG. 1, the x-ray tomographic system 100 that can be used in the 
present invention uses a gantry 110. The gantry 110 contains an x-ray point source 105 
that projects a beam 120 at a detector array 1 15 on the opposite side of the gantry 110. 

20 The beam 120 passes through the subject 125, and the individual detectors 130 within 
the detector array 115 sense the attenuation of the beam 120 passing through the subject 
125. The detectors 130 generate electrical signals corresponding to the attenuation, and 
the x-ray source and detector assembly rotates about the subject 120 to generate 
projection data. 

25 Detector array 1 15 is formed by detector elements or cells 130 which together 

sense the projected x-rays that pass through an object 125, for example a patient. 
Detector array 115 may be fabricated in a single slice or multi-slice configuration, and 
also, in various shapes, such as an arc of a circle, an arc of a cylinder, or flat panels. 
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The detectors 130 can transmit the projection data to the computer system 135 or other 
system processor or processors, which can reconstruct an image from the projection 
data. The computer 135 can, in some embodiments, transmit this image to a sub- 
system for display, and possible archiving, which in one embodiment can be the 
computer's display 135D. In one embodiment, the computer 135 can be a SiUcon 
Crraphics 02 computing platform (Silicon Graphics, Inc.; Mountain View, CA, USA), 
or any other suitable single or multiprocessor computing system such as those further 
described herein below. 

The image reconstruction systems and methods can be executed on one or more 
system processors. Additionally, the systems and methods can be programmed in any 
programming language that one skilled in the art would find appropriate. In one 
embodiment, the system can be programmed in IDL (Interactive Data Language from 
Research Systems, Boulder, CO, USA). 

Rotation of gantry 1 10 and the operation of x-ray source 105 are governed by a 
control mechanism of CT system 100. The control mechanism can reside inside or 
outside of computer 135 and can include an x-ray controller that provides power and 
timing signals to x-ray source 105 and any gantry motor controller that controls the 
rotational speed and position of gantry 110. A data acquisition system in the control 
mechanism can sample analog data from detector elements 130 and convert the data to 
digital signals for subsequent processing. In some embodiments, these digital signals 
can be stored in a system data store (SDS) described in further detail below. An image 
reconstructor using the me±ods described in attachments 1 and 2 can receive sampled 
and digitized x-ray data from the data acquisition system and perform image 
reconstruction. The reconstructed image can be applied as an input to a system data 
store or system processor, as further described hereinbelow, or computer 135. 

In one embodiment, computer 135 can also receive commands and scanning 
parameters from an operator via console that has a keyboard. An associated cathode 
ray tube display 135D, or other suitable display device (e.g., liquid crystal display, 
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plasma display, etc.), can be used to allow the operator to observe the reconstructed 
image and other data from computer 135. The operator supplied commands and 
parameters are used by computer 135 to provide control signals and information to data 
acquisition system, x-ray controller and gantry motor controller. 

According to an exemplary embodiment, the image processing and 
reconstruction system includes a system processor potentially including multiple 
processing elements. The term processing element may refer to (1) a process running 
on a particular piece, or across particular pieces, of processing hardware, (2) a 
particular piece of processing hardware, or either (1) or (2) as the context allows. Each 
processing element can be supported via a standard general purpose processor such as 
the processor within the Silicon Graphics 02 computing platform discussed 
hereinabove; alternative processors such as UltraSPARC (Sun Microsystems, Palo 
Alto, CA) or an Intel-compatible processor platforms preferably using at least one 
PENTIUM in, PENTIUM IV or CELERON (Intel Corp., Santa Clara, CA) class 
processor could be used in other embodiments. The system processor can include one 
or more field programmable gate arrays (FPGAs) and/or application specific integrated 
circuits (ASICs) configured to perform at least a portion of the functionality according 
to the present invention. In some embodiments, the system processor can include a 
combination of general purpose processors, ASICs and/or FPGAs. In some 
embodiments, image processing and reconstruction, as further described in attachments 
1 and 2, can be distributed across multiple processing elements. In some such 
embodiments, aspects of the functionality, or portions thereof, may be executed in 
series or in parallel; particular functionality, or portions thereof, executed a multiplicity 
of times may also occur in series or parallel. In FIG. 1, computer 135, or portions 
thereof, can be included within, or serve as, the system processor of the present 
invention. 

In a system processor including at least one general purpose processor, the 
general purpose processor typically runs an appropriate operating system such as 
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WINDOWS/NT, WINDOWS 2000 or WINDOWS/XP (Microsoft, Redmond, WA), 
Solaris (Sun Microsystems, Palo Alto, CA), or LINUX (or other UNIX variant). In 
some embodiments that include the Silicon Graphics 02 platform discussed 
hereinabove, the operating system used is Silicon Graphics' UNIX variant. 
5 The SDS referenced above may include a variety of primary and secondary 

storage elements. In one embodiment, the SDS may include RAM as part of the 
primary storage; the amount of RAM might range, for example, from 256 MB to 2 GB 
in embodiments including a computer workstation. One skilled in the art will recognize 
that depending on the image resolution, the amount of storage used can vary 
10 significantly. The primary storage can in some embodiments include other forms of 
memory such as cache memory, registers, non-volatile memory (e.g., FLASH, ROM, 
EPROM, etc.), etc. 

The SDS can also include secondary storage including single, multiple and/or 
varied servers and storage elements. For example, the SDS may use internal storage 

15 devices connected to the system processor. In embodiments where a single processing 
element supports all of the system functionality, such as computer 135 in FIG. 1, a local 
hard disk drive can serve as the secondary storage of the SDS, and a disk operating 
system executing on such a single processing element can act as a data server receiving 
and servicing data requests. A system bus would serve as the communication channel 

20 between the system processor and the SDS (typically, at least RAM and the hard disk 
drive). 

It will be understood by those skilled in the art that the different information 
used in the analysis and control processes and systems according to the present 
invention can be logically or physically segregated within a single device serving as 
25 secondary storage for the SDS; multiple related data stores accessible through a unified 
management system, which together serve as the SDS; or multiple independent data 
stores individually accessible through disparate management systems, which may in 
some embodiments be collectively viewed as the SDS. The various storage elements 
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that comprise the physical architecture of the SDS may be centrally located, or 
distributed across a variety of diverse locations. 

The architecture of the secondary storage of the system data store may vary 
significantly in different embodiments. In several embodiments, database(s) are used 
5 to store and manipulate the data; in some such embodiments, one or more relational 
database management systems, such as DB2 (IBM, White Plains, NY), SQL Server 
(Microsoft, Redmond, WA), ACCESS (Microsoft, Redmond, WA), ORACLE 8i 
(Oracle Corp., Redwood Shores, CA), Ingres (Computer Associates, Islandia, NY), 
MySQL (MySQL AB, Sweden) or Adaptive Server Enterprise (Sybase Inc., 

10 Emeryville, CA), may be used in connection with a variety of storage devices/file 

servers that may include one or more standard magnetic and/or optical disk drives using 
any appropriate interface including, without limitation, DDE and SCSI. In some 
embodiments, a tape library such as Exabyte X80 (Exabyte Corporation, Boulder, CO), 
a storage attached network (SAN) solution such as available from (EMC, Inc., 

15 Hopkinton, MA), a network attached storage (NAS) solution such as a NetApp Filer 
740 (Network Appliances, Sunnyvale, CA), or combinations thereof may be used. In 
other embodiments, the data store may use database systems with other architectures 
such as object-oriented, spatial, object-relational or hierarchical. 

Instead of, or in addition to, those organization approaches discussed above, 

20 certain embodiments may use other storage implementations such as hash tables or flat 
files or combinations of such architectures. Such alternative approaches may use data 
servers other than database management systems such as a hash table look-up server, 
procedure and/or process and/or a flat file retrieval server, procedure and/or process. 
Further, the SDS may use a combination of any of such approaches in organizing its 

25 secondary storage architecture. 

The SDS communicates with the system processor by one or more 
communication channels. Multiple channels can be involved in some embodiments for 
supporting conmiunication between processing elements of the system processor and 
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portions of the SDS. Such channels can include without limitation computer network, 
direct dial-up connection, dedicated connection, direct or indirect connection such as 
via a bus connection, parallel or serial connection, USB connection, null modem 
connection or wireless connection utilizing an appropriate communication protocol 
such as BLUETOOTH, IRDA, 802.1 lb or other suitable channel as would be known to 
those skilled in the art. 

All forms of data, including raw, intermediate, and computed can be stored on 
one or more SDS either temporarily or permanently. 

In one embodiment, a Grangeat-type half-scan algorithm is used in the circular 
scanning case. In another embodiment, the Grangeat-type half-scan algorithm is 
extended to work in a helical case without data truncation. In yet another embodiment, 
the Grangeat-type half-scan algorithm is extended to work in a helical case with data 
truncation. These embodiments are described in further detail below. 

A Grangeat-Type Half-Scan Algorithm for Cone-Beam CT 

According to a first exemplary embodiment, a Grangeat-type half-scan 
algorithm is provided in the circular scanning case. Appropriate data filling is 
performed using the Grangeat approach, and the off-mid-plane artifacts associated 
algorithms, such as the Feldkamp-type algorithms, are suppressed. 

In the following, Grangeat' s formula for circular half-scan geometry is modified 
by combining redundant data to improve contrast resolution while optimizing temporal 
resolution. The redundant data are weighted with smooth functions to suppress data 
inconsistence. 

As shown in FIG. 2, the Radon transform of a 3D-function fix) is defined by 




(1) 
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where n is the unit vector which passes through the characteristic point C described 
by spherical coordinates (p, 0,<p), x Cartesian coordinates (x, z) . Eq. (1) means 

that the Radon value at C is the integral of the object function /(3c) on the plane 

through C and normal to the vector n . It is well known that the 3-D function f{x) 

can be reconstructed from Rf{pn) provided that Rf(/On) is available for all planes 

through a neighborhood of point X . The inversion formula of the 3D Radon transform 



For cone-beam CT, it is instrumental to connect cone-beam data to 3D Radon 
data. One connection is Grangeat's formulation, which is geometrically attractive. 
Mathematically, as shown in FIG. 3, the link can be expressed as follows: 



is given by 




(2) 



^Rf{pn)^Ry{pn) 




(3) 




210064 



ATTORNEY DOCKET NO. 21087.0025U2 UIRF 02077 

PATENT 

16 

where Xf[s(fin)j^l/r{/On)] is the detector value the distance ^ away from the detector 

center O along the line t perpendicular to on the detector plane located at 

the angle ^ from the y-axis, denotes the distance between the source and the 

origin, the distance between the source and an arbitrary point A along ^ , and ^ 
5 denotes the angle between the line SO and SC, and 

SO 

X^f[s(pn),t,vripn)] = — Xf[s(pn\t,t^^ 

Given a characteristic point ^ in the Radon domain, the plane orthogonal to the 
vector pn is determined. Then, the intersection point(s) of the plane with the source 

trajectory ^^^^ can be found, and the detector plane(s) specified, on which the line 

C D 
10 integration can be performed. Let ^ denote the intersection of the detector plane ^ 

with the ray that comes from ^^"^^ and goes through ^ , The position can be 
described by a vector sHp . To compute the derivative of Radon value at ^ , the line 
integration is performed along ^ , which is orthogonal to the vector sfi^ . 

For a digital implementation of Eq. (3), the derivative in the s -direction is 
15 reformulated as the sum of its horizontal and vertical components as Eq. (4) below: 



d 3 d 

— X ^f[si/Dn), t, yf{fjn)\ = cos[a(/»i)]— X^f[sipn),t, wifi^)} + sm[a(fm)]— X ^f[s{fin\t, yrifln)] 
as dp dq 

where p and q are the Cartesian axes, s and a define a polar system on the detector 
20 plane (See FIG. 3(c)). Substituting Eq. (4) into Eq. (3), we have Eq. (5) as follows: 

Rjifm) = — \—{cos[a{fM)] X ^fUiph), t,\ir{fm)\dt + sin[a(pi)] X^f[s{piiX t, y/{fmy\dt 

a SO 
where cos p = 
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Grangeat-Type Half-Scan Formula 

Eq. (3) can be modified into the following half-scan version. 

d ^ 13" SO 

— Rf{pn) = Y,Q^iifin) \— Xf(s(pn),t,l//,(pn))dt (6) 

O/) COS OS ijA 

5 where 

Vr, (p, 0,(p) = (p-^ sin (7) 
L50 sm^J 

^, «:^) = ^ + ^ - sin "M J^, (8) 
L5Osin0J 

and a)^{p,0,(p) and a}^ip^B^(p) are smooth weighting functions explained in detail 
10 below. 

A derivation of the Grangeat-type half-scan algorithm follows. 

Derivation of the Grangeat-type half-scan algorithm 

As indicated above, Grangeat found a link between Radon data and cone-beam 
15 projection data, which is expressed as follows. 

^Rf(pn) = R'fipn) = — ]^Xf[s(pnUMpn)}dt 
dp cos p ds _iSA 

(al) 

= -^-^3- \x^f[sipnU,Wipnm 
cos p as _i 

In the circular full-scan case, there always exist two detector planes from which the 
Radon data can be calculated except for the shadow zone. The detector planes are 
specified by the following two equations: 

20 vr^ (p, 0,^) = (p-^ sin — ^— , (a2) 

ISO sin 0 J 
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In an ideal situation where data is noise-free and the object is stationary during 
the rotation, we have Eq. (a4) which follows: 

5 

Therefore, the radon data can be calculated from either of the detector planes. 
However, in practice we should consider the following practical conditions: (1) 
Projection data contain noise; (2) Data from D^^ and D^^ may be different due to 

motion of an object. 
10 Hence, we have: 

Considering noise in projection data, we should combine redundant data to improve the 
signal-to-noise ratio. Consequently, Eq. (a5) can be modified to the following: 

■^Rfipn) = — W|- ^\^Xf{s{pn),t,ii/,{pn))dt, (a6) 

uf) i-\ COS p OS _^ o/l 

15 where (O^ (fin) + CO2 (fin) = 1 . To maximize the signal to noise ratio, both fi^i (fin) and 

^2(/^) should equal to 1/2. 

However, when motion is significant, one of the solutions to suppress the 

motion artifacts is to use half-scan, which is an important goal according to exemplary 

embodiments. But in this case, data are not always redundant by a factor of two. More 
20 exactly, there are generally doubly, singly sampled regions, and shadow zones on a 

meridian plane. Therefore, there exist discontinuities between the adjacent regions. 

Hence, the smooth weighting functions such as those described below are needed to 

suppress the associated artifacts. 
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Exemplary half-scan geometry is depicted in FIG. 4, which is the central plane 
where the source trajectory resides. The cone angle is denoted as , defined as a half 

of the full cone-angle. The scanning angle ^ varies from 0 to ;r + 2y^ . The 
horizontal axis of a detector plane is denoted as p . 

In the circular full-scan case, for any characteristic point not in the shadow zone 
or on its surface there exist a pair of detector planes specified by Eqs. (7) and (8). 
However, in the circular half-scan case, the dual planes are not always available. When 
the dual planes are found, we are in a "doubly sampled zone"; when one of them is 
missing due to the half-scan, we are in a "singly sampled zone". Of course, there is not 
an associated detector plane in the shadow zone. Relative to the full-scan, the shadow 
zone is increased due to the half-scan. The shadow zone area on a meridian plane 
depends on the angle (p of the meridian plane. 



Boundaries between Singly and Doubly Sampled Zones 

The boundary equations between these zones are instrumental for design of 
weighing functions Q)i(p,0,<p) and Q}^i^p^O^(p) , as well as interpolation and 

the angular difference between the meridian plane M ^ and the detector plane , 

where the line integration is performed. This angular difference ranges from -90 to 
+90 degrees. Given a half-scan range and a meridian angle (p , it can be determined 



extrapolation into the shadow zone. The term sin ^ 



sin^J 



in Eqs. (7) and (8) is 



whether D^^ and D^^ are available by evaluating sin M- 



It is emphasized 



that the following geometric property is important to understand the relationship 
between singly and doubly sampled zones. For a given characteristic point, the normal 
directions of the two associated detector planes must be symmetric to the normal 
direction of the meridian plane containing that characteristic point. 
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Denoting A{(p) as the critical angle which separates the singly and doubly 
sampled zones, as shown in FIGS. 5(a)-(d), then the number of the available detector 
planes changes from two to one or one to two across this critical angle. Specifically, 
for the meridian plane of an angle q> the boundary between the singly and doubly 



5 sampled zones is expressed as sin . 

' "Osin^ 



-if Pb 
[50 sir 



= A{g>) . Therefore, we have 



(0, (p) = 50 sin A{(p) sin 6 . (9) 
The critical angle function A{q)) takes different forms depending on the 
meridian angle <p . We have the following four intervals: 

U:0<ip<r^, (10) 
10 12: r„<<p<90, 

13: 90<^<90 + 2;^^, 
14: 90 ^2y^ <<p<180. 

If a pair of detector planes, D^j^^^^ and D^+^^yi(^) , are available, there exist two 
redundant data for {p,6,q>) , If only either of them is available, there exists only one 
15 data. If none of them is available, there exists no data. In FIGS. 5(a)-(d), the 

availability of D^^^^^^^ and D^^^.^^^^ are marked showing whether both are available, 
only either of them available, or none of them available for a given meridian plane 
M ^ . There are the four intervals in FIGS. 5(a)-(d). For II, both ^^+^(^) and 

are available when -(p< A((p) <nl2. Only D,^^^^^^^^ is available when 
20 2y„ -q>< A{q)) < -q> , Therefore, the boundary between doubly and singly region is 
formed at A;i (jp) = —<p . Hence, boundary equation becomes 

(0, (p) = SO sin(-^) sin 0 from Eq. 9. In the same manner, the critical angle 
functions can be acquired for the other intervals and are expressed as follows: 

II: A,,((p) = -(p, (11) 
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12: A,^{(p)^<p-lYm^ 

13: A,3(^) = ^-2;^„,andA,3.(^) = «?-:^, 

14: A^^{(p) = (p-n . 

These critical angle functions are plotted in FIG. 6. It should be noted that there 
5 exist two critical functions for 13. Therefore, there exist two boundaries in this interval, 
given as {0, <p) = SO sin( (^)) sin 6 and p[ (0, <p) = SO sin( A^y (^)) sin 0 , 
respectively. 



Shadow Zone Boundaries 

10 The shadow zone boundaries are needed to interpolate/extrapolate for missing 

data. The boundary equation for the shadow zone in the full-scan case is given by 
=50sin«9) (12) 
In other words, each detector plane gives a Radon circle of diameter SO on the 
meridian plane. However, in the half-scan case, the diameter of the Radon circle may 
15 change, depending on the meridian plane angle, as shown in FIGS. 7(a)-(b). 

Specifically, given a meridian plane, we can find both the diameters of the two Radon 
circles on the plane. Then, we can express the shadow zone boundaries as follows: 

= SO, m sin(^) , p0<O, (13) 
p,=SO^(<p)sm(0), p0>O. 
20 In reference to FIGS. 7(a)-(b), it is easy to verify that 

Xo=SO/2, =0, (14) 
X, =(50/2)cos(2^„+;r), 
=(50/2)sin(2r„+;r). 
Then, it can be geometrically derived in each of the four intervals that 

25 II : SO, «p) = N ^^'^\ ^ so, ((p) = SO ; (15) 

cos ecan 
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12: so,,<p)J^^^lf^^^,SOA<P)-SO; 
[COS ec<p\ 

13: 50; {(p) = SO , SO, {<p) = 50 ; 

\2xy -2>, cot^l 



UIRF 02077 



14: SOi((p) = SO, SO^((p) = ^ 



cos ecc 



Weighting Functions 

S Using the above boundary equations, we can graphically depict the singly and 

doubly sampled zones, along with the shadow zones. Four maps from different 
meridian planes are represented in FIGS. 8(a)-(d). White area stands for the doubly 
sampled region, gray for the singly sampled region, and black for the shadow zone. 

When data are consistent, we can set (O^—Q)^^^ for the doubly sampled region, set 

10 the valid one to 1 and the other as 0 for the singly sampled region, and set both to zero 
for the shadow zone where missing data will be estimated afterwards. When data are 
inconsistent in practice, the following smooth weighting functions are designed 
according to Parker's half-scan weighting scheme: 
II: 



15 (O^{p,0,(p) = cos 



^[n p-{-S) 

0, 
0, 



sm' 



oi^{^p,0,(p) - sin 



^2S-p,my 



,U p-i-5) ' 

yi p^{e)-{-S)j 



1, 
1, 



9<0 , p<p,i0); 

0<O , p>p,(0); 
0>O, p<p,i0); 

a>0, p>p,(0). 

0<O, p<p^i0); 

0<O , p>p,(0); 
0>O, p<p^i0); 



(16) 



(17) 
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cos' 



10 



15 



7t p-p^iei " 
2 s-p,{e)j 



12: 



cOi(p,0,(p) = sin 



cos 



(o^ip,d,(p) = cos 



sin 



13: 



cOjip,0,<p) = sin' 



cos 



sin 



^ p-{-S) \ 
p,{d)-{,-5)) 

1, 
1, 

'tc p-p,{ef 
^2S-p,(.e)J' 

J Jt p-{-5) 

yi p.m-i-s)^ 

0, 
0, 

Jn_p-p^ 

s-p.my 

'n p-{-S) ^ 
^2p,(0)-(-S)^ 

1, 

;r p-Pi, (0) 
^2s-p,\0)) 

r 

It p- 
1, 



^>0, p>p,(0). 



0<(i, p<p,{0); (18) 

^<o , p>p,{d)- 

0>O, p<p,i0); 
0>O, p>p,i0). 

0<O,p<p,{0); (19) 

0<O , p>p,(0); 
0>O, p<p,i0); 

0>O, p>p,i0). 



0<O, p<p^(0); (20) 



0<O , p,i0)<p<p, {0); 



0<O, p>p, (0); 



0>O, p<p, (J9); 



0>O, p, (0)<p<p^{0}; 
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O)^{p,0,0) = COS' 



sin 



cos' 



sm' 



7t p-{-S) 

{2p,m-i-s)) 

0, 



n P-Ph W) 

{'^s-p^e)] 

n p-{-5) 
0, 

'n p-p,{ei ^ 

,2 5-p,(e)y 



14: 

(0^{p,e,q>) = \, 



cos 



sin' 



co^{p,e,q>)=0. 



sin 



cos 



n p-p^{ef 
2 S-p,{e)^ 

p-i-S) ^ 
2 P,(0)-i-S) 

1, 



'n p- plies ' 
y2S-p,{e)y 

Jtt p-i-S) \ 

\2p,{e)-{-S)y 



24 



e>Q, p>p,{9). 

e<Q, p<p^{e); (21) 

d<o , p^{e)<p<p^{d)- 
e<o, p>pl{e); 

e>o, p<p^{e)\ 

6>>0, p,\0)<p<p,i0); 

e>o, p>p,{e). 

e<o,p<p,m; (22) 
d<o, p>p^ie); 

e>0, p<p,(0); 

0>O, p>p„{0). 
0<O,p<p,m\ (23) 

0<O, p>p,i0); 
0>O, p<p^m\ 
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0, ^ 0>O,p>p,(0). 

In the above weighting functions, S is assunaed to be the half size of the p 
support and ps [-S,S] . Some representative weighting functions are shown in FIGS. 9 
(a)-(b). 

5 

Interpolation/Extrapolation 

According to exemplary embodiments, zero-padding and linear interpolation 

methods may be used, respectively, to estimate missing data in the shadow zone. 

Grangeat reconstruction after zero padding is theoretically equivalent to the Feldkamp- 
10 type reconstruction. A few of Radon values near the shadow zone boundary can be 

linearly interpolated along 0 direction. The zero padding and linear interpolation are 

shown in FIGS. 10(a) and 10(b). 

Of course, there are multiple possibilities for interpolation/extrapolation into the 

shadow zone. Knowledge based interpolation/extrapolation is also feasible in this half- 
15 scan Grangeat framework. In addition, parallel-beam approximation of cone-beam 

projection data can be applied in this rebinning framework as well. 

Exemplary Implementation of the Grangeat-type half-scan algorithm for the 
circular scanning case 
20 To summarize, the Grangeat-type half-scan algorithm can, according to a first 

embodiment, be implemented in the following steps: 

(1) Specify a characteristic point (p,^,^) where the derivative of Radon data can be 
calculated; 

(2) Determine the line integration points for the given characteristic point according to 
25 Eqs. (7), (8), and appropriate rebinning equations; 

(3) Calculate the derivatives of Radon data; 

(4) Calculate the smooth weighting functions a)^{p,0^(p) and (O^^P-'^^V) \ 
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(5) Apply the weighting functions to the Radon data; 

(6) Repeat Steps 1-5 until we are done with all the characteristic points that can be 
calculated; 

(7) Estimate the Radon data in the shadow zone; and 

5 (8) Reconstruct an image according to a Radon inversion formula. 

Exemplary Simulation Results 

Using a software simulator in the BDL Language (Research Systems Inc., 
Boulder, Colorado) for Grangeat-type image reconstruction, the numerical 

10 differentiation was performed with a built-in function based on 3-point Lagrangian 
interpolation. The source-to-origin distance was set to 3.92. The number of detectors 
per cone-beam projection was 256 by 256. The size of the 2D detector plane was 2.1 
by 2.1. A half of the full-cone angle is about 15 degree. The number of projections 
was 210. The number of meridian planes was 180. The numbers of radial and angular 

15 samples were 256 and 360, respectively. Each reconstructed image volume had 
dimensions of 2.1 by 2.1 by 2.1, and contained 256 by 256 by 256 voxels. The 
numbers of samples were chosen to be greater than the lower bounds we established 
using the Fourier analysis method. Both the spherical phantom and the 3D Shepp- 
Logan phantom as shown in Table 1 were used in the numerical simulation. 

20 Table 1 . Parameters of the phantoms used in our numerical simulation. 



Phantom 


a 


b 


C 


X 


y 


z 


e 


<p 


Density 


Sphere 


0.9 


0.9 


0.9 


0.0 


0.0 


0.0 


0.0 


0.0 


2.0 


Shepp- 
Logan 


0.69 


0.9 


0.92 


0.0 


0.0 


0.0 


0.0 


0.0 


2.0 




0.6624 


0.88 


0.874 


0.0 


0.0 


-0.0184 


0.0 


0.0 


-0.98 




0.41 


0.21 


0.16 


-0.22 


-0.25 


0.0 


0.0 


72.0 


-0.02 




0.31 


0.22 


0.11 


0.22 


-0.25 


0.0 


0.0 


-72.0 


-0.02 




0.21 


0.35 


0.25 


0.0 


-0.25 


0.35 


0.0 


0.0 


0.01 
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0.046 


0.046 


0.046 


0.0 


-0.25 


0.1 


0.0 


0.0 


0.01 




0.046 


0.02 


0.023 


-0.08 


-0.25 


-0.605 


0.0 


0.0 


0.01 




0.046 


0.02 


0.023 


0.06 


-0.25 


-0.605 


0.0 


90.0 


0.01 




0.056 


0.1 


0.04 


0.06 


0.625 


-0.105 


0.0 


0.0 


0.02 




0.056 


0.1 


0.056 


0.0 


0.0625 


0.1 


0.0 


0.0 


-0.02 




0.046 


0.046 


0.046 


0.0 


-0.25 


-0.1 


0.0 


0.0 


0.01 




0.023 


0.023 


0.023 


0.0 


-0.25 


-0.605 


0.0 


0.0 


0.01 



FIG. 1 1 shows the results obtained from the y=0 plane of the sphere phantom. 
The half-scan algorithm cannot produce exact reconstruction for a general object due to 
the incompleteness of projection data but in special cases like a sphere phantom this 
5 algorithm can achieve exact reconstruction with linear interpolation. When data in the 
shadow zone are linearly interpolated from the half-scan data of a sphere, exact 
reconstruction of the sphere can be achieved because the radon transform of a perfect 
sphere is linear. This exactness seems impossible with other half-scan cone-beam 
reconstruction algorithms. Nevertheless, this property is desirable in a number of major 
10 applications, such as micro-CT studies of spherical specimens. 

FIG. 12 presents the results obtained from the y=0.242 and x=-0.0369 planes of 
the 3D Shepp-Logan phantom. With the zero-padding method, the low intensity drop 
away from the center plane was found to be more serious than that with the Feldkamp 
half-scan reconstruction. With the linear interpolation method, this type of artifacts 
15 was substantially eliminated. 

FIG. 13 compares the full-scan and half-scan Grangeat algorithms. Linear 
interpolation was applied with both algorithms. As compared to the Grangeat 
reconstruction in the full-scan case, little image degradation was visually perceived in 
all the Grangeat-type half-scan reconstructions performed. 
20 FIG. 14 compares Feldkamp-type half-scan reconstruction and Grangeat-type 

half-scan reconstruction. 
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In the first exemplary embodiment described above, a Grangeat-type half-scan 
algorithm is used in the circular scanning case. The half-scan spans 180 degrees plus 
two cone angles, allowing exact in-plane reconstruction. The smooth half-scan 
weighting functions have been designed for suppression of data inconsistence. 
5 Numerical simulation results have suggested the correctness and demonstrated the 
merits of our formulas. This Grangeat-type half-scan algorithm is considered 
promising for quantitative and dynamic biomedical applications of CT and micro-CT. 



Grangeat-Type Hettcal Half-Scan CT Algorithm for Reconstruction of a Short 
10 Object 

According to a second embodiment, image reconstruction may be performed in 
the helical case without data truncation. In this embodiment, Grangeat's formula may 
be modified for utilization and estimation of Radon data. Specifically, each 
characteristic point may be characterized in the Radon space into singly, doubly, triply 

15 sampled, and shadow regions, respectively. A smooth weighting strategy is designed to 
compensate for data for redundancy and inconsistency. In the helical half-scan case, 
the concepts of projected trajectories and transition points on meridian planes are 
introduced to guide the design of weighting functions. Then, the shadow region is 
recovered via linear interpolation after smooth weighting. The Shepp-Logan phantom 

20 and other phantoms can be used to verify the correctness of the formulation, and 

demonstrate the merits of the Grangeat-type half-scan algorithm. The Grangeat-type 
helical half-scan algorithm is not only valuable for quantitative and/or dynamic 
biomedical applications of CT and micro-CT, but also serves as a basis in solving the 
long object problem. 

25 

Rebinning Equations for a Helical Trajectory 

The geometry for a helical half-scan is shown in FIG. 15, where (x,y,z) is the 
reference coordinate system, {y\z') denotes a meridian plane at ip from the x-axis, and 
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the source vertex, A. , ranges from 0 to n-^-ly^ . If the helical pitch is denoted by /i , the 
source trajectory can be parameterized as 



To compute the derivative of the Radon value at >aw , we need first calculate the 
line integration point C£> on a detector plane, through which a line integration is 
performed along the line normal to OCp . There is a geometrical relationship between 
the characteristic point (/7,^,^)and the line integration point (j. a, A) . In other words, 
we can find {s,a,A) from (p,^,^). 

The 3D Radon value at a characteristic point ip,d,<p) is the integration of an 
object on a plane satisfying the following equation: 

hmx^p. (25) 
The intersection of the plane with the source trajectory is found by replacing x with 
a{A) to solve 

nma(A) = p (26) 

where n = (sin ^ cos ^, sin ^ sin ^, cos ^) denotes a unit vector towards the corresponding 
characteristic point, and a(A) = {R cos /I. R sin X, hA-^zo) the source position. The equation 
can be written as: 



While we can calculate the line integration point analytically in a circular trajectory 
case, we can only do it numerically in a helical trajectory case. Once we acquire A , we 

have line integration point (a , 5 ) on the detector plane , where i// = A-\-^. The 

equations are expressed as: 



diA) = (SO cos A, SO sin A,hA + Zo)- 



(24) 



R cos A sin 9 cos (p+ RsinA sin ^ sin ^ + (hA + Zq ) ^^s 6- p = Q 
Rsin ^(cos/lcos^ + sin yisin^) + (/i/l + Zo)cos^ = p 



(27) 
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£ir = tan (■ 



cot^ 



(28) 



sin(^->l) 



s = 



(p - {hA + zq ) cos 0) sin a 
cosO 



(29) 



Therefore, the Radon value at the characteristic point defined by (p, 0, cp) can be 
calculated by the integration along the line represented by (5, a, A) according to Eqs. 
27-29. 

Grangeat-Type Helical Half-Scan Formula 

With a circular half-scan, there are three types of regions: shadow, singly and 
doubly sampled regions, respectively. With a helical half -scan, in addition to those 
three types of regions, there may be triply sampled regions as well. They are 
schematically illustrated in FIGS. 18(a)-(d). Hence, the Grangeat-type helical half-scan 
formula must be in the following format: 



where (o^ denote the weighting functions, ^ is a group identifier to be explained in 
detail later, and 



The form is basically the same as that with a circular half-scan trajectory but three 
weighting functions and corresponding Radon values are needed for each characteristic 
point, while we only need two weighting functions in a circular half-scan case. 

In Eq. (31), the value of the group identifier g is determined according to the 
following criteria: 
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g = 



1, 



0<A<A 



(32) 



3, 



/l,2 <A<^-\-2y^ 



where A^^ and tangential vertices and will be explained in detail below. This 

means that the calculated Radon data must belong to one of the three groups depending 
on the A . Regarding the weighting functions, they must satisfy 



If we assume that there is no motion or data inconsistency during the scan, we 
can simply average Radon data from different vertices. In this case, the weighting 
functions become discontinuous. However, in real situations, motion effects and data 
inconsistency should be taken into account. Then, the discontinuous weighting 
functions could cause artifacts. Therefore, the smooth weighting functions are needed 
for satisfactory image quality, as explained below. 

Weighting Functions 

The concepts of projected trajectory and transition points are needed before our 
weighting functions are designed. 

Projected Trajectory on a Meridian Plane at (p 

The projected trajectory on a meridian plane at ^ is useful for design of the 
weighting functions. Several projected trajectories on different meridian planes are 
shown in FIGS. 19(a)-19(i). The source trajectory of Eq. (24) is first transformed from 



3 




(33) 
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(jc,y,z) to (x\y\z') coordinate systems, as shown in FIG. 17. Therefore, for a given 
meridian plane at ^, the projected trajectory on this plane is described as 



Dashed lines in FIG. 20(a) represent the planes normal to the meridian plane. The 
integration result(s) on the plane(s) corresponding to each line must be equal to the 
Radon value at (p,0,<p) . In the Grangeat framework, the derivative of Radon data is 
calculated from the detector planes as denoted by the colored lines associated with the 
intersected vertex points. Geometrically, there can be maximally three intersection 
points. For a given {0, <p) , the number of intersecting points determines the degree of 
redundancy and depends on p . 

Transition Points on the Projected Trajectory 

Every projected trajectory has two end points expressed by 



For a given ^ on a meridian plane <p , there may be planes normal to 3 and 
tangent to the projected trajectory. The tangential points can be analytically specified 
by 




(34) 




(35) 



(36) 




(37) 



(38) 
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where An , A,^ are calculated in the following. 

A projected trajectory on a meridian plane is determined as 

y = SOsina-{(p~)), (39) 



To express A as a function of z , we obtain 



(40) 



h 



(41) 



substitute Eq. (43) into Eq. (41) , and have 

y = SOsin(^^-(<p-^)). (42) 
n 2 

Then, we compute the derivative and set it to the slope of the colored line: 

10 = -50cos(^-^-(^--)) = -cot^ . (43) 

dz h h 2 

In other words, 

z = Kcos^' (- cot ^) + - f )) + zo . (44) 
Then, we have 

X = (cos-k-^cot^) + (^-|)) . (45) 
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The solution of this equation is meaningful only when 0<A.<fr+2y„ , resulting 
in up to two solutions: A^i and A,2 ■ It is assumed that A,2 is greater than A,i if it exists. 
Recall that these and A,2 are used for grouping vertices as described above. 

In terms of the above end points and tangential points, we can find the transition 
5 points p as follows: 

Pel =^»^. 

(46) 

Pn =^•'1 
Pt2=^*h 

10 where 6 = (>, z) = (sin 0, cosB) . As shown in FIG. 20(a) the type of a region can be 
identified according to the following criterion: p^i<p< p^i , indicating a doubly 
sampled region; Pe2<P< Pe\ » a singly sampled region; p^^<p< p,, , a doubly sampled 
region. 

IS Smooth Weighting Functions 

Weighting functions are designed in reference to p^i , p^2 > Pa > Pti » which are 
functions of 9 and cp . Mathematically speaking, there are up to two possibilities when 
we have neither p^, nor Pf2 , which are p^^ ^ Pe2 and p^^ - Pei • Similarly, there are up to 
six cases when we have only either >9,i or p,2 • Furthermore, there are up to twenty-four 
20 cases when we have both and Pt2 • However, all the statements are not meaningful 
after our case-by-case inspection. It is found that there exist only eleven cases actually. 
For example, in absence of p^i and » we always have the case of p^i < p^2 > and never 
have the case of p^i > p^2 • Similarly, it is impossible to have the cases of p^, < Pfi < p^2 
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and p^2 - Pt\ - Pel absence of • Those eleven cases are: i) p^x - Pei absence of 
Ptx and Pfi ; ii) Pt\ ^ Pe\ ^ Pe^ absence Pa \ iii) Pel - Pei - Pa absence of \ iv) 
< p^2 ^ Pel in absence of p,^ ; v) /?^, < p^2 ^ Pa in absence of p,2 \ vi) 

Pt2 ^ Pel ^ Pel ^ Pa ; vii) Pn ^ Pel ^ Pel ^ Pt2 \ viu) Pel ^ Pti ^ Pa ^ Pel \ ix) Pt2 < Pel ^ Pel ^ Pa ; 

X) /?,2<Pei<Ai<P,2;xi) P,^<P,2<Pe2^Pn- 

In addition to our geometric analysis on projected trajectories on the meridian 
plane, we did exclusive numerical simulation to confirm that there are only the above 
eleven meaningful cases. For all the iO,<p), we calculated the transition points, sorted 
the p values, and decided which one of the mathematically possible thirty two cases it 
belongs to. Once a specific case was found, we set the flag for that case on. After this 
kind of numerical verification, we eliminated the cases that never happened. Also, we 
repeated the simulation with respect to representative combinations of imaging 
parameters, including the source-to-origin distance, helical pitch and cone angle ( 2y^ ). 
Finally, it was confirmed that there are indeed only the above eleven cases that are 
feasible. 

For each of those meaningful cases, smooth weighing functions are designed in 
terms of p^i , p^2 > Pa > Pa • Th^ two general requirements are that (1) the sum of the 
weighs for each characteristic point be one, and (2) the weight profile along the p 
direction be continuous. To further understand the designing processing, some 
representative illustrations are considered helpful. FIG. 20(b) represents the case where 
there is neither p^^ nor Pa for given (^,^), and only group 1 exists. Therefore, =1 , 
and and co^ are set to zero. FIG. 20(c) corresponds to the case where there is one 
tangential point, and two groups are available. The trajectory from to F, belongs to 
group 1, while that from F, to belongs to group 2. The weighting function for 
groupl, co^ , is designed to change gradually from 1 to 0 along the projected trajectory 
from p„ to , while the weighting function for group 2, (O2 , increases from 0 to 1 in 
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the same interval, and stays constant between p^, and p^i • ^3 is set to zero since there 
is no group 3. Of course, the sum of the weighting functions should be made one. 
FIG. 20(d) illustrates the case where there are two tangential points, and we have group 
1 between and F, , group 2 between \ and , and group 3 between F^ and ^2 ■ 
5 weighting function for the group 1, , should be one between p^x to A2 > and smoothly 
change from 1 to 0 along the projected trajectory from to p^i . The weighting 
function for the group 3, co^ , should smoothly change from 0 to 1 along the trajectory 
from /9,2 to Pfx^ and be 0 between p^^ and p^2 • Th^ weighting function for the group 2, 
(O2 , is designed to smoothly increase from 0 at p^i until it reaches the middle point 
10 between p^ and p,i , and decrease to 0 at p^i . 

Some representative distributions of the weighting functions are included in 
FIGS. 21(a)-(l). The weighting functions are formulated as follows, keyed to each of 
the feasible cases: 

i) < absence of p^x and p^^ • 
15 a)i =1, 

Q)^ =0. 

ii) Pa < Pel ^ Pel in absence p^ : 

cos^C^-^::^) , Pa^P<Pei 

^ Pel ~Pt\ 

20 0 , p,x<P<Pel 

^ Pel ~Pt\ 

1 , p,x<P<Pe2 

co^= 0 
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iii) /?e2 - Pa - Pa absence of /?,2 '• 

<w, = 0' , p^2^P<Pa 

sin^(ff^) , PA^P^Pa 

^2= 1 . Pe\^P<Pe2 

5 cOS^(^^Z^) , P.,<P<Pa 

0)^= 0 

iv) iO^i < p^2 - Pe\ absence of • 

^ sin2(£^Z^) , Pa^P^Pel 

1 , p^2<P^Pe\ 

10 ^2= COS^(~ ^"^^^ ) , Pa^P^Pel 

2 Pel- Pa 

0 , p^2<P^Pe\ 

0)^= 0 

v) < ^ Ai in absence of p,2 • 

2 Ai - A2 



0)2= 0 , p^^<p<P^2 

"2 Pa- Pel 



sin2(£^_^) , Pei^P^Pn 



0)^= 0 

Vi) p,2<Pe2^Pel^Pil 
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a>i= 0 
0 



2 Ai -P,i 



0)2 = sin^(- ^ ^'^ ) 

2 P«2 - A2 



1 



COS ( ^- ^-^) 

2 Ai - Ai 



2,^ P-Pn ^ 
O)^ = cos^(- '-^) 

2 Pe2-P,2 



0 
0 

Vii) p,i<p^^<p^2^P,2 

2 p„ -Ai 

0 
0 



15 sin^(^^^::^) 

2 Prf -P,i 



2 P,2-Pt2 



a>3= 0 



38 

P,2<P<P,2 
Pe2<P<Pel 

Pa^P^P,\ 

P,2<P<P^2 

Pe2 <P<Pa 
P,l<P<P,^ 

P,2<p<p^2 

Pe2 <P<Pa 
PA^P^PtX 



Pn^P^Pel 

Pel <P<P.2 
Pe2^P^Pi2 

Pn^P^PA 

Pa <P<Pe2 

P^2^P<P,2 

Pix^P^Pa 
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15 



0 

Sin ( — - — 

2 Pa- Pel 

viii) p,i<p,2<p,i<p,2 

0}y= 1 



cos^(;zr ) 

Pl\ -Pt2 



0 
0 

0)2= 0 



-:-2^„ P-P,2 s 

Sin ) 

Pn -Pti 



Ai -Pa 



0 

<Mj= 0 

0 



Pt\ -Pa 



1 

ix) p,2<P,x<P^2<p,i 

<y, = 0 



isi„^(£-£::^) 

2 2 



39 

/Orf <P<P,2 

p,2<P^Pa 
Pei^P<Pa 

Pa^P<Pa+iPii-Pa)l'2- 
Pa + ^P,i-Pa)l'^^P^Pn 

Pt\<P^Pe2 

P.\^P<Pa 

p,2<p<Pa+(fi,i-Pa')l'^ 
Pa-^iPti- Pa)l'^^ P^ P,\ 

P,l<P^Pe2 

Pe\^P<Pa 

p,2<p<Pa+iP,\-Pa)l'2' 
Pa+iPii- Pa)l'^^ P^ Pt\ 

P,\<P^Pe2 

P,2<P<Pa 
Pa ^P^Pei 
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1 

^ 2 

2 

j_ 
2 



^2= 2 



0 



) 



fi;, = 0 



Pt\-Pe\ 



0 
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Pel<P^Pt\ 
Pn^P<Pe\ 
Pel ^P^ Pel 

P,l^P<Pe\ 

Pel ^P^Pel 
Pe2<P^P,l 

Pa<P<pA 

Pel^P<PelHPa-pel)l^ 

Ai+(Ai-Ai)/2<p<Ai 

Pn<P^Pe2 
P,2^P<Pel 

A.i^P<P*i + (Ai-Ai)/2 
Ai+(Ai-Ai)/2</?<Ai 

P,l<P^Pe2 

P,2<P<Pa 
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15 



icos^(;r^^) 
2 Ai -Pa 

cos^(;z- ) 
P,\ -Pa 

1 

xi) p,^<Pa<P,2<P,x 

a)i= 1 

cos^(;r 

Pe2-Pi2 
2 Pe2-P,2 

2 

a>2 = 0 
0 



icos^(;r-^:^) 

2 Pe2 -P,2 



(0,= 0 



Pel -P,2 



Pe2 -Pi2 



0 



41 

Pa^P<PaHP,i-Pa)I'2- 
/'rf+(Ai-Prf)/2<p</>„ 

P.l<P^Pe2 

Pa^P<P,2 

P,2<P< P,2 + iPe2 -Pi-dl"^ 

P,2 + iPe2 - Aa) / 2 < iO < p^2 

p,2<P<P,i 
PA^P<Pt2 

Pa^P< P,2 + (.Pe2 - Pr2) / 2 
P,2 + iPe2-Pi2)''^^P^Pe2 

P,2<P<P,i 

Pe\^P<Pa 

P,2<p<p,2+(.P,2-P,2)l'2 
Pt2 + (Pe2 - >£?,2) / 2 < P < 
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Although examples of weighting functions have been described above, it should 
be appreciated that the invention is not limited to these weighting functions. Other 
weighting functions may also be used, designed in substantially the same spirit of the 
design of the weighting functions described above. 

5 Interpolation 

As in the circular scanning case described above, linear interpolation may also 
be used in the helical scanning case to estimate missing data. The boundaries of the 
shadow zone were numerically determined. Then, the shadow zones were linearly 
interpolated along the d direction using the measured boundary values. To 
10 demonstrate this interpolation idea graphically, FIGS. 22(a)-(b) show the derivatives of 
Radon data with no interpolation and with linear interpolation, respectively. 

It should be appreciated that other interpolation methods may also be used, 
designed in substantially the same spirit of the interpolation procedure described above. 

15 Exemplary Implementation Procedure 

To summarize, the Grangeat-type half-scan algorithm with a helical scanning 
locus can be implemented in the following exemplary steps: 

(1) Specify a characteristic point 9^ (p) where the derivative of Radon data can be 
calculated; 
20 (2) Calculate X,^ and \ 

(3) Calculate Pa^Pti.Pei^Pei '^ 

(4) Calculate smooth weighting functions a)^(p,e,<p) , co^ip.O.q)) , and^yaCp.^,^) ; 

(5) Determine line integration points for the given characteristic point according to 
appropriate rebinning equations; 

25 (6) Calculate the derivatives of Radon data and store them according to their group 
membership as determined, according to one embodiment, by Eq. (23); 
(7) Apply the weighting functions to the derivative of Radon data; 
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(8) Repeat Steps (l)-(7) until all the measurable characteristic points are done; 

(9) Estimate the Radon data in the shadow zone (using the linear interpolation method 
according to one embodiment); and 

(10) Use an inverse Radon formula to reconstruct an image volume. 

5 According to an exemplary embodiment, results from steps (2)-(4) can be pre- 

calculated and stored for computational efficiency. Similarly, the rebinning 
coefficients from Step (5) can be found beforehand. 



Exemplary Results 

10 A software simulator was developed in the DDL Language (Research Systems 

Inc., Boulder, Colorado) for Grangeat-type image reconstruction. In the 
implementation of the Grangeat-type formula, the numerical differentiation was 
performed with a built-in function based on 3-point Lagrangian interpolation. 

The source-to-origin distance was set to 5. The number of detectors per cone- 

15 beam projection was 256 by 256. The size of the 2D detector plane was 3.3 by 3.3. 

The heUcal pitch was 2. The full-cone angle was about 30 degree. The scan range was 
from 0 to ^ + 2y„ . The number of projections was 210, The number of meridian planes 
was 180. The numbers of radial and angular samples were 256 and 360, respectively. 
Each reconstructed image volume had dimensions of 4.1 by 4.1 by 4.1, and contained 

20 256 by 256 by 256 voxels. 

One might use the same number of projections above and below each and every 
slice so that all slices have the same image quality. In one approach, an asymmetric 
number of slices above and below the slice were used, except for the z=0 slice. If we 
use the symmetric half-scan helical Feldkamp, which has symmetric projections for any 

25 z location, all slices would have similar image quality. However, it means that each 
slice is reconstructed in a different time window. 

According to this embodiment, a half-scan algorithm is provided in the 
Grangeat framework with superior temporal characteristics (temporal resolution and 
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temporal consistency; the latter requires that a whole volume of interest is reconstructed 
with projections in the same time window) and less image artifacts (which is achieved 
by appropriate data handling in the shadow zone). The reason that the asymmetric 
projections were used was to maintain this temporal consistency. Half-scan helical 
5 Feldkamp methods may also be used in the same way. 

The 3D Shepp-Lx)gan phantom was used in the numerical simulation as shown 
in Table 2. 



Table 2. Parameters of the phantoms used in our numerical simulation. 



Phantom 


a 


B 


C 


X 


y 


z 


e 


(p 


Density 


Shepp- 
Logan 


0.69 


0.9 


0.92 


0.0 


0.0 


0.0 


0.0 


0.0 


2.0 




0.6624 


0.88 


0.874 


0.0 


0.0 


-0.0184 


0.0 


0.0 


-0.98 




0.41 


0.21 


0.16 


-0.22 


-0.25 


0.0 


0.0 


72.0 


-0.02 




0.31 


0.22 


0.11 


0.22 


-0.25 


0.0 


0.0 


-72.0 


-0.02 




0.21 


0.35 


0.25 


0.0 


-0.25 


0.35 


0.0 


0.0 


0.01 




0.046 


0.046 


0.046 


0.0 


-0.25 


0.1 


0.0 


0.0 


0.01 




0.046 


0.02 


0.023 


-0.08 


-0.25 


-0.605 


0.0 


0.0 


0.01 




0.046 


0.02 


0.023 


0.06 


-0.25 


-0.605 


0.0 


90.0 


0.01 




0.056 


0.1 


0.04 


0.06 


0.625 


-0.105 


0.0 


0.0 


0.02 




0.056 


0.1 


0.056 


0.0 


0.0625 


0.1 


0.0 


0.0 


-0.02 




0.046 


0.046 


0.046 


0.0 


-0.25 


-0.1 


0.0 


0.0 


0.01 




0.023 


0.023 


0.023 


0.0 


-0.25 


-0.605 


0.0 


0.0 


0.01 



10 In Table 11, parameters of the phantoms used in our numerical simulation, (a, b, 

c) denote the semi-axes of ellipsoids, (x, y, z) the center coordinates of each ellipsoid, 
6 and <p the same as in FIG 15. The actual density at a position is determined by 
sununing the densities of the ellipsoids covering that point. 



210064 



ATTORNEY DOCKET NO. 21087.0025U2 
PATENT 



UIRF 02077 



45 

FIG. 23 shows the derivatives of Radon data of the Shepp-Logan phantom, the 
weighting functions for each group, and the combined data. FIGS. 24(a)-(c) present 
typical reconstructed sHces of the Shepp-Logan phantom. With the helical half-scan 
Feldkamp method (FIG. 24(a)) and the helical half-scan Grangeat and zero-padding 
5 method (FIG. 24(b)), the low intensity drop was serious away from the center plane. 
However, this type of artifacts was essentially eliminated with the helical half-scan 
Grangeat and linear interpolation method (FIG. 24(c)). 

According to an exemplary embodiment, parallel-beam approximation of cone- 
beam projection data (in the spatial domain) may also be used to estimate missing data. 

10 In addition to the spatial domain approximation, the Radon domain estimation, such as 
linear interpolation, spline interpolation and knowledge-based interpolation, can be 
done in our framework as well. However, in the filtered backprojection framework, 
each frame of cone-beam projection data can be processed as soon as it is acquired, a 
desirable property for practical implementation. 

15 According to exemplary embodiment, the strengths of the Radon space based 

algorithm described above may be combined with an algorithm for achieving high 
computational efficiency, e.g., that reported by Noo and Heuscher in the Deftise-Clack 
framework, described in Noo, F. and D.J. Heuscher, Image Reconstruction from cone- 
beam data on a circular short-scan, SPEE Medical Imaging, 2002, herein incorporated 

20 by reference. In this algorithm, a parallel-beam approximation is used to estimate 
missing cone-beam data due to the incompleteness of the half-scan. 

The Radon space based algorithm yields more flexible data filling and can be 
modified for long object reconstruction, while the filtered backprojection algorithm 
provides high computational efficiency. Logically, we expect that in development of 

25 the Radon space based algorithm we will find the optimal strategy for missing data 
recovery; then we will seek for a filtered backprojection implementation of that 
algorithm without compromising the best achievable image quality. Our proposed 
Altered backprojection algorithm is similar to what has been described in the preceding 
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subsection, except that we will apply a Noo-Heuscher-type filtered backprojection 
method for the short object reconstruction after using one of the following two filtered 
backprojection algorithms for the first order reconstruction: 

The helical scanning geometry described herein solves the short object problem, 
which assumes that the object is completely covered by the X-ray cone beam from any 
source position. Even though research with the short object geometry is valuable on its 
own, such as for micro-CT imaging of spherical samples, the geometry for data 
truncation along the axial direction is the most practical. It should be appreciated that 
the Grangeat-type helical half-scan CT work may also be extended to the long object 
case, as described below. 

According to yet another embodiment, the Grangeat-type half-scan algorithm 
described above may be used for reconstruction of a long object, which has important 
applications for clinical CT and biomedical micro-CT. 

According to this embodiment, a temporally non-optimized pre-reconstruction 
step may be used to transform the long object problem into a short object problem. The 
detector area is analytically classifled into desirable, corrupted, and useless areas. The 
cone-beam data in the corrupted area is corrected by the forward projection of the pre- 
reconstruction data, while data in the useless area is set to zero. Wang's generalized 
Feldkamp-type algorithm and Schaller's PHI algorithm may be used for the pre- 
reconstruction. 

The Wang generalized Feldkamp algorithm is described, e.g., in Wang, G., et 
al., A General Cone-beam reconstruction algorithm, IEEE Trans. Med. Imaging, 1993, 
12: pp. 483-496, herein incorporated by reference. This algorithm is approximate but 
efficient. Its half-scan version allows high temporal resolution but suffers from shading 
artifacts for a large cone angle. 

The Schaller PHI algorithm is described, e.g., in Schaller, S., et al., Exact radon 
rebinning algorithms for the long object problem in helical cone-beam CT, IEEE Trans. 
Med. Imaging, 2000. 19(5), pp. 361-375. Tam proposed to scan along a helix and two 
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circles, and avoid the truncated data in the rebinning framework. The Schaller PHI 
algorithm does not need the troublesome circles. In the Schaller method, Radon data 
on each meridian plane are calculated from the virtual object bounded by the helix lines 
instead of the two circles. Note that Radon data from the PHI algorithm are temporally 
poor due to the fact that the data are synthesized from several helical turns. 

After the correction, the half-scan Grangeat method for the short object is 
applied along with several shadow zone interpolation techniques to correct the image 
artifacts. Specially designed temporal resolution phantoms are adopted to assess the 
image quality over the volume of interest. Comparison with a half-scan Feldkamp 
algorithm is conducted. 

The Wang algorithm is the first effective yet simplest attempt to the long object 
problem, but the nature of the algorithm is only approximate and can be considered as 
the worst estimation. In our simulation, we achieved very similar image quality as we 
do in the short object case even with the worst estimation for the pre-reconstruction. 
The Schaller*s method is very compatible with our general scheme and facilitates the 
Radon domain shadow zone interpolation. In conclusion, our proposed Grangeat-type 
half-scan strategy for solving the long object problem produces clearly superior image 
quality than Feldkamp-type algorithms, although its temporal resolution and image 
artifacts are substantially affected by the cone angle. 

Exemplary Implementation Procedure 

In general, the steps involved for long object image reconstructions are as 
follows (See FIG. 25). 

(1) Select one of existing long object algorithms (none of which is optimized for the 
image quality we need, especially in terms of temporal resolution); 

(2) Specify a volume of interest, which is the short object to be reconstructed; 

(3) Perform the first-order reconstruction using the algorithm selected in (1). 
According to one embodiment, there are two possibilities: (3.1) reconstruct relevant 
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regions outside the short object, and (3.2) reconstruct relevant regions inside the short 
object; 

(4) Classify measured cone-beam data into (4.1) desirable, (4.2) corrupted, and (4.3) 
useless data (Fig. 7(a)). The desirable data come exclusively from the short object, 
while the corrupted data involve structures inside and outside the short object, and the 
useless data carry no information on the short object; 

(5) Correct the corrupted data based on forward projection of relevant voxel values 
in the first-order reconstruction. There are the two ways corresponding to the two 
possibilities in (3): (5.1) subtract segment integrals outside the short object from the 
measured cone-beam data, and (5.2) directly use segment integrals inside the short 
object; 

(6) Estimate missing data on the shadow zone in the Radon space using forward 
projection of the first-order reconstruction and/or interpolation of Radon data 
derivatives calculated from measured cone-beam data; 

(7) Perform necessary featuring between the data from the direct measures and the 
data from the. first-order reconstruction; 

(8) Perform our Grangeat-type half-scan reconstruction to reconstruct the short 
object for high temporal resolution. 

The modification of our proposed algorithm into the filtered backprojection 
format will greatly improve its computational complexity. 

Exemplary Simulation Results 

Simulation of this embodiment revealed that image quality of the long object 
reconstruction is similar to that of the short object reconstruction. The Schaller's 
method is very compatible with this scheme, gives the benefits in the Radon domain 
data manipulation, and produces slightly better image quality than that with Wang's 
algorithm. The image artifacts are affected by the cone angle but are much less serious 
than that with the Feldkamp-type algorithm. 
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Throughout this application and attachments hereto, various publications may 
have been referenced. The disclosures of these publications in their entireties are 
hereby incorporated by reference into this application in order to more fully describe the 
state of the art to which this invention pertains. The publications disclosed are also 
individually and specifically incorporated by reference herein for the material contained 
in them that is discussed in the sentence in which the reference is relied upon. 

The embodiments described above are given as illustrative exaniples only. It 
will be readily appreciated by those skilled in the art that many deviations may be made 
from the specific embodiments disclosed in this specification without departing from 
the scope of the invention. For example, although the embodiments described above 
have been directed to circular and helical scanning cases, it should be appreciated that 
the invention is also applicable to substantially circular and substantially helical, spiral, 
and spiral-like scanning cases. 
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